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We simulate head-on collisions from rest at large separation of binary white dwarf - neutron stars 
(WDNSs) in full general relativity. Our study serves as a prelude to our analysis of the circular 
binary WDNS problem. We focus on compact binaries whose total mass exceeds the maximum mass 
that a cold degenerate star can support, and our goal is to determine the fate of such systems. A fully 
general relativistic hydrodynamic computation of a realistic WDNS head-on collision is prohibitive 
due to the large range of dynamical time scales and length scales involved. For this reason, we 
construct an equation of state (EOS) which captures the main physical features of NSs while, at 
the same time, scales down the size of WDs. We call these scaled-down WD models "pseudo-WDs 
(pWDs)". Using these pWDs, we can study these systems via a sequence of simulations where the 
size of the pWD gradually increases toward the realistic case. We perform two sets of simulations; 
One set studies the effects of the NS mass on the final outcome, when the pWD is kept fixed. The 
other set studies the effect of the pWD compaction on the final outcome, when the pWD mass and 
the NS are kept fixed. All simulations show that after the collision, 14%- 18% of the initial total 
rest mass escapes to infinity. All remnant masses still exceed the maximum rest mass that our cold 
EOS can support (I.92M0), but no case leads to prompt collapse to a black hole. This outcome 
arises because the final configurations are hot. All cases settle into spherical, quasiequilibrium 
configurations consisting of a cold NS core surrounded by a hot mantle, resembling Thorne-Zytkow 
objects. Extrapolating our results to realistic WD compactions, we predict that the likely outcome of 
a head-on collision of a realistic, massive WDNS system will be the formation of a quasiequilibrium 
Thorne-Zytkow- like object. 

PACS numbers: 04.25.D-,04.25.dk,04.40.Dg 



I. INTRODUCTION 

The inspiral and merger of compact binaries represent 
some of the most promising sources of gravitational waves 
(GWs) for detection by ground-based laser interferome- 
ters like LIGO [l|,[l, VIRGO Bi, GEO [5], TAMA 
and AIGO as well as by the proposed space-based 
interferometers LISA and DECIGO [10|. Extracting 
physical information from gravitational radiation emitted 
by compact binaries and determining their ultimate fate 
requires careful modeling of these systems in full general 
relativity (see [11] and references therein). Most effort 
to date has focused on modeling black hole-black hole 
(BHBH) binaries (see [12] and references therein), and 
neutron star-neutron star (NSNS) binaries (see [13] for 
a review), with some recent general relativistic work on 
black hole-neutron star (BHNS) binaries p!3 - [32| . 

In this work we consider white dwarf-neutron star 
(WDNS) binaries in full general relativity. WDNS bi- 
naries are promising sources of low-frequency GWs for 
LISA and DECIGO and, as we argued in [33], possi- 
bly also high-frequency GWs for LIGO, VIRGO, GEO, 
TAMA and AIGO, if the remnant ultimately collapses to 
form a black hole. 

Like NSNS binaries, WDNS binaries are known to ex- 
ist. In [33] we compiled tables with 20 observed WDNS 
binaries and their orbital properties. The NS masses 
in these systems range between I.26M0 and 2.O8M0, 
and their distribution is centered around l.SM©. On 



the other hand, the WD masses in these systems range 
between O.125M0 and I.3M0, and their distribution is 
centered around O.GM©. Finally, 18 of these observed 
WDNS binaries have total mass greater than l.GSM©, of 
which 8 have a WD component with mass greater than 
O.8M0, and 5 have total mass greater than 2.2M0. This 
is interesting because the expected TOV limiting mass 
for a cold, degenerate gas ranges between l.GSM© and 
2.2M0 [34-42], depending on the equation of state (EOS) 
and degree of rotation, and one of the main goals of this 
work is to determine whether a WDNS merger can lead 
to prompt collapse to a black hole. 

Population synthesis calculations by Nelemans et al. 
m show that there are about 2.2 x 10^ WDNS bina- 
ries in our Galaxy, and that they have a merger rate of 
1.4 X 10~^yr~^. Furthermore, Nelemans et al. find that 
after a year of integration, LISA should be able to de- 
tect 128 WDNS pre-merger binaries and, after consider- 
ing the contribution of the double WD background GW 
noise, resolve 38 of these. On the other hand, calcula- 
tions by Cooray [H] give much more conservative num- 
bers of resolvable WDNS binaries. In particular, Cooray 
finds that the number of LIS A-resolvable WDNS binaries 
ranges between 1-10, using a WDNS merger rate between 
10~^yr~^-10~^yr~^. Finally, recent work by Thompson 
et al. [45] suggests that the lower limit on the merger rate 
in the Milky Way, at 95% confidence, is 2.5 x 10~^yr~^. 
Thompson et al. also suggest that the merger rate in the 
local universe is ~ 0.5 — 1 x 10^Gpc~^yr~^. Therefore, 
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leaving some uncertainties aside, all recent work on pop- 
ulation synthesis suggests that LISA should be able to 
detect a few WDNS pre-mergers per year. 

We note here that Newtonian work on binaries with a 
WD component has been performed analytically in [33> , 
|46[45QI1^ and via Newtonian hydrodynamic simulations in 

In [33I we focused on WDNS binaries that have spiraled 
sufficiently close that they reach the termination point 
for equilibrium configurations. This is the Roche limit 
for WDNSs, at which point the WD fills its Roche lobe 
and may experience one of at least two possible fates: 
i) stable mass transfer (SMT) from the WD across the 
inner Lagrange point onto the NS, or ii) tidal disruption 
of the WD by the NS via unstable mass transfer (UMT). 

We also studied the key parameters that determine 
whether a system will undergo SMT or UMT and found 
that, for a given NS mass, there exists a critical mass 
ratio g'crit ^ 2/3 which separates the UMT and SMT 
regimes. If the mass ratio q = Mwd/^ns of a WDNS 
system is such that q > ^crit, the WD quickly overfills 
its Roche lobe, and the binary will ultimately undergo 
UMT. In the opposite case, q < ^crit, the system will un- 
dergo SMT. We showed that a quasistationary treatment 
is adequate to follow the evolution of an SMT binary dur- 
ing this secular phase and calculated the gravitational 
waveforms. We also pointed out that WDNS observa- 
tions suggest that there are candidates residing in both 
regimes. 

In the case of tidal disruption (UMT), by contrast, the 
system will evolve on a hydrodynamical (orbital) time 
scale. In this scenario the NS may plunge into the WD 
and spiral into the center of the star, forming a quasiequi- 
librium configuration that resembles a Thorne-Zytkow 
object (TZO) [57]; alternatively, the NS may be the re- 
ceptacle of massive debris from the disrupted WD. 

Depending on the details of the EOS, a cold degenerate 
gas can support a maximum gravitational mass between 
1.65M0 and 2.2M0 against catastrophic collapse, if it is 
not rotating (the TOV limit). It can support 20% more 
mass, if it is rotating uniformly at the mass-shedding 
limit (a "supramassive NS" [36]), and about 50% more 
mass, if it rotates differentially (a "hypermassive NS" 
j33-l36|). If the total mass of the merged WDNS exceeds 
the maximum mass supportable by a cold EOS, delayed 
collapse to a black hole is inevitable after the remnant 
cools. However, the ultimate fate of the merged WDNS 
depends on the initial mass of the cold progenitor stars, 
the degree of mass and angular momentum loss during 
the WD disruption and binary merger phases, the an- 
gular momentum profile of the WDNS remnant and the 
extent to which the disrupted debris is heated by shocks 
as it settles onto the NS and forms an extended, massive 
mantle. These are issues that require a hydrodynamic 
simulation to resolve. Moreover, ascertaining whether or 
not the neutron star ultimately undergoes a catastrophic 
collapse (either prompt or delayed) to a black hole re- 
quires that such a simulation be performed in full gen- 



eral relativity. In fact, even the final fate of the NS in 
the alternative scenario in which there is a long epoch 
of SMT may also lead to catastrophic collapse, if the 
neutron star mass is close to the neutron star maximum 
mass, and this scenario too will require a general rela- 
tivistic hydrodynamic simulation to track. 

In this paper we employ the Illinois adaptive mesh re- 
finement (AMR) relativistic hydrodynamics code [23,^581 
to perform our first simulations of these alternative sce- 
narios. In particular, we study the head-on collision from 
rest at large separation of a massive WD and a NS as a 
prelude to our investigation of the circular binary prob- 
lem, which we will report in a future work. We focus on 
compact objects whose total mass exceeds the maximum 
mass supportable by a cold EOS to determine whether 
such a collision leads to prompt collapse of the remnant, 
or a hot gaseous mantle composed of WD debris sur- 
rounding a central NS - a Thorne-Zytkow-like object 
(TZIO). 

The vast range of dynamical time scales and length 
scales involved in this problem make fully general rela- 
tivistic simulations extremely challenging. For example, 
a near-Chandrasekhar-mass WD has a radius Rwb — 
lO^km and dynamical time scale of about Is. On the 
other hand a typical NS has a radius of order R^s — 
10km and dynamical time scale of about 1ms. Therefore, 
there is a difference of several orders of magnitude both 
in length scales and time scales. Current numerical rel- 
ativity techniques and available computational resources 
make such calculations prohibitive. For this reason, we 
tackle this problem using a different strategy. 

In particular, we construct a piecewise polytropic EOS 
which captures the main physical features of a NS while, 
at the same time, scales down the size of the WD. We 
call these scaled-down WDs "pseudo-WDs (pWDs)". We 
perform a sequence of simulations where we change the 
EOS so that the pWDs have the same mass (O.98M0) but 
different compactions, while the compaction and mass of 
the NS involved remain practically unchanged. In other 
words, while keeping the masses of the binary compo- 
nents and the NS radius fixed, we adjust the ratio of the 
radius of the pWD to that of the NS so that it varies from 
5:1 to 20:1 and then use our results to predict the out- 
come of the realistic case. The common feature among all 
versions of the piecewise EOS we employ is that the max- 
imum NS mass always is I.8M0 and the maximum WD 
mass always is 1.43M0, i.e., the Chandrasekhar mass. 

In addition to studying the effects of the pWD com- 
paction, we also study the effects of the NS mass. We 
consider NSs with masses 1.4M0 I.5M0 and I.6M0. 

All simulations that we perform show that after the 
collision, 14%- 18% of the initial total rest mass escapes to 
infinity. The remnant mass in all cases exceeds the maxi- 
mum rest mass that our cold EOS can support (I.92M0), 
but we find that no case leads to prompt collapse to a 
black hole. This outcome arises because the final con- 
figurations are hot. All our cases settle into a spher- 
ical quasiequilibrium configuration consisting of a cold 
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NS core surrounded by a hot mantle. Hence, all rem- 
nants are TZlOs. Extrapolating our results to realistic 
WD compactions, we predict that the likely outcome of 
a head-on collision from rest at large separation of a re- 
alistic massive WDNS system will be the formation of a 
quasiequilibrium TZIO. 

This paper is organized as follows. In Section |II] we 
review the time scales and length scales involved in a 
WDNS merger and discuss why this problem presents 
such a computational challenge. In Section |TTl] we in- 
troduce the EOS adopted for our computations and de- 
scribe our pWD models. Section HVl outlines our method 
for preparing initial data, and Section IVl summarizes our 
methods for evolving the gravitational and matter fields. 
We present the results of our fully relativistic hydrody- 
namic simulations in Section IVH and conclude in Sec- 
tion |VlI] with a summary of our main findings. Through- 
out we use geometrized units, where G = c = 1. 



II. COMPUTATIONAL CHALLENGE 

Simulating a WDNS merger in full general relativity is 
a difficult computational task. In this section we sketch 
in quantitative terms exactly why this is so. 

There are three fundamental time scales and length 
scales involved in the WDNS merger that must be re- 
solved. The relevant time scales are the dynamical time 
scale of each component of the binary and the orbital pe- 
riod; the relevant length scales are the NS and WD radii 
and their orbital separation. 

Resolving the WD length scale and dynamical time 
scale is necessary in order to assess what happens to the 
WD at merger. Merger occurs on the orbital time scale, 
so this time scale must also be resolved. Resolving the NS 
dynamical time scale will enable us to assess whether the 
NS promptly collapses and forms black hole, or remains 
inside the remnant WD, settling into a TZIO. 

The dynamical time scale of the NS, td,NS5 is given by 



^d,NS 



Mns'' 



(1) 



where i^NS^ and M^s are the characteristic NS radius 
and mass respectively. 

Similarly, the dynamical time scale of the WD, td,wD, 
is given by 



^d,WD 



WD 



(2) 



where Rwb , and Mwd are the characteristic WD radius 
and mass respectively. Finally, the orbital time scale, T, 
is given by 



T = 27ri 



Mr' 



(3) 



where Mt = Mns + Mwd is the total mass, and A is 
the orbital separation. Note that, at this separation, the 
head-on collision time scale is 



^coll 



A3 



2^2 V^T 4^2' 



(4) 



assuming the stars free-fall from rest as Newtonian point 
masses, and hence it is roughly the same order of magni- 
tude as T. 

By use of Eqs. ([T]) and ^ the ratio of the WD time 
scale to the NS time scale is 



^d,WD 
^d,NS 



qR^ 



■NS 



where 



Mwd 
Mns 



(5) 



(6) 



is the binary mass ratio. 

For Mwd = IM© and using a cold degenerate electron 
EOS one finds i^wD ^ 5000km [H, El. On the other 
hand, typical NS masses and radii are Mns = l.SM©, 
^NS ^ 10km. Hence, in realistic scenarios the ratio of 
the WD size to the NS size is 



^NS 



500, 



(7) 



and from Eq. (|5]) the ratio of the dynamical time scales 
is 



^™ ~ 10^ 

^d,NS 



(8) 



At the Roche limit, A is typically a few (two to five) 
times the WD radius [33]. Using, A ~ 2i?wD, and the 
values for the masses and radii used above, the orbital 
time scale becomes 



T 

^d,NS 



A^ 



(1 



q)Rf 



10^ 



(9) 



NS 



It is thus clear that there is a vast range of length scales 
and time scales involved in this problem. The only way 
to simulate the WDNS merger is by exploiting the power 
of adaptive mesh refinement, so that resolution is high 
only where required. However, even this does not suffice 
to tackle the timescale problem as we explain below. 

Given that all current numerical relativity schemes for 
evolving both the spacetime and the ffuid are explicit, 
there are strong limitations imposed on the size of the 
timestep by the Courant-Friedrich-Lewy condition 



At_ 

Ax 



= A<C, 



(10) 



where A is the Courant number and C a constant of order 
unity that depends on the integration scheme employed. 
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If one uses AMR, this implies that the size of the timestep 
has to be different in regions of different mesh size Ax. 
If we resolve the stars adequately, the mesh size will be 
much smaller near the NS than near the WD, because 
in typical scenarios the NS is 500 times smaller than the 
WD. Eq. (fTQ|) then implies that the smallest timestep 
must be in the domain of the NS. In particular, if the 
NS is covered by A/'ns = 2i?NS / Axns grid zones and the 
WD is covered by A^wd = 2i?wD/AxwD grid zones, then 
from Eq. (p!Q]) we have 



AtwB AxwD A/'ns ^wd 



At 



NS 



Axns Awd ^ns 



(11) 



where A^ns^^^ns and AtwD,AxwD denote the 
timestep and mesh size in the vicinity of the NS and 
WD, respectively. 

To assess the potential formation of a black hole re- 
quires at least A/'ns = 50 grid zones across the NS, in 
order that a 2Mq BH (which is a probable mass a BH 
would have after the merger of a WDNS system of to- 
tal mass of about 2.5Mq) would be covered by at least 
20 grid zones. Even if a BH does not form, covering 
the NS with 50 grid zones is necessary to reliably model 
the NS and maintain small hamiltonian and momentum 
constraint violations. Resolving the WD requires about 
A^WD = 30 grid zones to reliably model the star. If we 
combine Eq. ([7]) and Eq. (fTT]) . we obtain 



AtwD 
AtNS 



833. 



(12) 



This means that for one timestep in the vicinity of the 
WD we would have to take about 833 timesteps in the 
vicinity of the NS. Even evolving the system for only 
one WD dynamical time scale would require millions of 
timesteps in the vicinity of the NS. This shows how diffi- 
cult it is to resolve both the WD and the NS at the same 
time. 

However, what renders the computation of the WDNS 
merger in full GR prohibitive is that a realistic merger 
takes place on an orbital time scale, which is equivalent 
to 10^ NS dynamical time scales (see Eq. (j9])). 

To make this quantitative, let us compare the orbital 
time scale with a typical timestep in the vicinity of the 
NS. Using A^NS grid zones across the NS and combining 
Eq. (dD and Eq. 1^ we find 



^d,NS A/'ns / ^NS 



AtNS 



2A V Mns 



100, 



(13) 



where in the last step we used a typical value for 
the Courant number A = 0.4 and the values for 
^ns^^ns^A'ns we cited above. Combining Eq. ([9j) and 
Eq. (p!3|) we obtain 



A^NS 



10^ 



(14) 



Hence, a realistic WDNS simulation would require a min- 
imum of 10'' timesteps in the vicinity of the NS. In fact 
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FIG. 1: Mass (M) - central rest-mass density (po,c) relation- 
ship of single TOV stars for various cold EOSs. In the plot 
Pnuc = 2 X 10^^ g/cm^ is the nuclear density. Plotted are the 
Chandrasekhar electron-degenerate EOS for mean molecular 
weight per electron /ie = 2 (labeled as WD), the AP2 version 
of the Akmal-Pandharipande-Ravenhall EOS [37, 60], a poly- 
tropic approximation of these realistic EOSs using EOS (|15p 
(labeled as Fit) and a version of EOS (p^ where the ratio of 
the isotropic radius of a 0.98Mq pWD to the isotropic radius 
of a 1.5M0 NS is reduced to ^ 10 (labeled as 10:1). The 
parameters of these EOSs are listed in Table [H 



this number of timesteps is an underestimate because 
extracting GWs would require a few orbits and the final 
system would settle in equilibrium or collapse within a 
few orbital time scales after merger. As a result, a dy- 
namical, fully general relativistic hydrodynamics WDNS 
calculation would require of order 10^ timesteps. 

We can give an estimate of how long such a simula- 
tion would be based on high resolution (192^ grid points 
in the innermost refinement level) benchmark runs we 
performed for a WDNS system with a I.SMq-NS and 
a 1.OM0-WD at separation of about 2.7 Rwb (close to 
the Roche limit), which has an orbital period of about 
1.4 X lO^M, where M = 2.5M0. Using 256 cores on the 
Ranger cluster of the Texas Advanced Computing Cen- 
ter we found that the Illinois GR hydrodynamics code 
advances about 6M per hour. Thus, the entire simula- 
tion (of about 10 orbital periods) would require about 
264 years of pure computational time. 

Realistic WDNS simulations are beyond the capabil- 
ities of current computational resources and numerical 
relativity techniques. For this reason, we will tackle the 
problem of WDNS mergers and head-on collisions adopt- 
ing an alternate strategy. We carefully construct an EOS 
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which mimics a reahstic cold NS EOS and, at the same 
time, scales down the size of the WD to make such a 
calculation feasible. Using sequences of these systems, 
where the WD size gradually increases, we can extrapo- 
late our results to the realistic case. These scaled down 
WDs or pseudo-WDs are the subject of the following sec- 
tion. 



IV. INITIAL DATA 

In this section we present the basic formalism for gener- 
ating valid general relativistic initial data for the head-on 
collision in a given pWDNS system. 

A. Gravitational Field Equations 



III. PSEUDO- WHITE DWARFS 

In this section we introduce our EOS and describe re- 
sulting models for pWDs. Our EOS is the following 6- 
parameter piecewise polytropic EOS 



P 

Po 



l/ni 



Po ^ Pi 



i^2pI^''\ Pi < Po ^ p2 



1/^3 

1^3 Po , P0> P2 



(15) 



where P is the pressure, po is the rest-mass density 
and /^i, /^2, /^3, ^1, ^2, ^3, Pi, P2 are the parameters of the 
EOS. Note that the parameters are 8 in number but con- 
tinuity requires that the following conditions be true 



1^1 



l/n2 — l/ni 
l^2pi 



1^2 = f^3p2 • (16) 

As a result, the adopted EOS ([15]) has 6 free parameters 
and throughout this paper we use a^s, ni, n2, na, pi, p2 to 
specify an EOS. Note also that we use the polytropic 
indices rii instead of the adiabatic indices Ti = 1 + l/n^. 

The freedom of our multi-parameter EOS enables us 
to capture the same characteristic curves and turning 
points on a TOV mass-central density plot as for a cold- 
degenerate realistic EOS (see [59]), as shown in Fig. [H 
The figure shows that EOS ([15]) can provide a rea- 
sonable approximation to the mass-central density rela- 
tion of realistic compact objects, exhibiting both stable 
{dM/dpo^c > 0) and unstable {dM/dpo^c < 0) branches 
for both WDs and NSs. 

Furthermore, EOS (fT5]) allows us to adjust the size of 
a pWD of given mass, thereby shifting the pWD branch 
to smaller radii (see Fig. [2]), while keeping the NS masses 
and radii approximately unchanged for Mns ^ 1.3Mq. 
The shifted branches in Fig. |2] correspond to the stars 
that we call pseudo-WDs in this work. 

Finally, note that all versions of EOS ([T5]) considered 
in this work have been carefully constructed so that the 
maximum gravitational mass of a NS is l.SM©, i.e., the 
same as that which for the AP2 version of the Akmal- 
Fandharipande-Ravenhah (APR) EOS [3?,'60], and the 
maximum gravitational mass of a pWD is 1.43M0, i.e., 
the maximum mass of a TOV WD obeying the Chan- 
drasekhar EOS for pe = 2. In addition, EOS (p!5]) is 
constructed to preserve the shape of the M vs po,c curve, 
yielding both stable and unstable branches and turning 
points appropriately. 



We begin by writing the spacetime metric in the stan- 
dard 3+1 form [61] 

ds^ = -a^dt^ + -fij{dx' + p'dt){dx^ + p^dt), (17) 

where a is the lapse function, the shift vector and 
the three-metric on spacelike hypersurfaces of constant 
time t. Throughout the paper Latin indices run from 1 
to 3, whereas Greek indices run from to 3. 

We conformally decompose the three-metric 7^^ as 



(18) 



where ^ is the conformal factor and fij the conformal 
metric. We adopt the standard approximation of a con- 
formally flat spacetime, so that fij = Sij in Cartesian 
coordinates. 

Since we are interested in head-on collisions between 
compact objects, we assume that initially the stars begin 
to accelerate towards each other starting from rest. As a 
result the extrinsic curvature is initially zero and the mo- 
mentum constraints are identically satisfied [62]. Hence, 
we need only prepare initial data for ^. 

Under the aforementioned assumptions the only equa- 
tion we have to solve is the Hamiltonian constraint, which 
becomes 

V^^ = -27r^V, (19) 

where is the flat Laplacian operator. Here the source 
term p is defined as 

p = n"n^T,^, (20) 

and where is the normal vector to a t = constant slice, 
and Tq,/3 is the stress-energy tensor of the matter. 

The gauge is chosen so that the initial slice is maximal, 
i.e., K = and dtK = 0, and the shift vector is set 
equal to zero. Using the assumption of maximal slicing 
it is straightforward to derive an equation for the lapse 
[11,^] 



V'a = 27ra^^{p^2S), 



where 



(21) 



(22) 



(23) 



and the source term S is defined as 

Equations ([T9|) and ([2T]) are elliptic and hence have to 
be supplemented with outer boundary conditions. Fol- 
lowing [63], we impose 1/r fall-off conditions on a — 1 
and ^ — 1 at the outer boundary. 
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TABLE I: Parameters for the piecewise poly tropic EOS ()15p used in generating different stellar models. The first column 
corresponds to the name of the EOS. An EOS named M : N corresponds to a version of ()15|) for which the mass - radius 
relationship of TOV stars is such that the ratio of the isotropic radius of a O.OSM© pWD to that of a l.SM© NS is M : N (see 
Fig. El). The EOS named AP2 is the same as the AP2 EOS defined in [6^. Finally, the EOS named Fit is an approximate fit 
to the Chandrasekhar EOS (for /ie = 2) joined onto the AP2 EOS (see Fig. [1]). Here pnuc = 1.48494 x 10~^km~^ and k,3 is 
given in geometrized units. 



EOS Name 


1^3 


ni 


712 


ns 


log(pi/pnuc) 


I0g(p2/Pnuc) 


20:1 


5064.2599 


1.56128 


2.98418 


0.714286 


-3.17219 


0.180473 


10:1 


4993.0688 


1.51515 


2.96971 


0.714286 


-2.26862 


0.208502 


5:1 


6123.5567 


1.51883 


2.94291 


0.699301 


-1.2909 


0.267623 


AP2 


145414.043 


0.60864 


0.49652 


0.514139 


0.398915 


0.698922 


Fit 


4458.0491 


2. 


2.96736 


0.716 


-6.39356 


0.208502 



B. Matter fields 



1. Diagnostics 



To model the matter, a perfect fluid stress-energy ten- 
sor is assumed: 



{po + Pi + P)u''u^ + Pg"P , 



a/3 



(24) 



where g^^ is the inverse of the four-metric and 
po^ pi^ P^u'^ are the rest-mass density, internal energy 
density, pressure, and four- velocity of the fluid respec- 
tively. 

Since the initial configuration is assumed to be at rest 
in the center of mass frame, the initial fluid four velocity 
is given by 



u"" = ^1^(1, 0,0,0) 



or 



(25) 



(26) 



A straightforward calculation shows that the source 
term p in Eq. (|2Q|) can then be written as 



and S in Eq. 



as 



P = Po ^ Pi, 



S = 3P. 



(27) 
(28) 



C. Computational methods 



To solve the elliptic equations (fT9|) and (|2T]) we de- 
veloped a fixed-mesh-refinement (FMR) finite difference 
code based on the Portable, Extensible Toolkit for Scien- 
tific Computation (PETSc) algorithms [64-66] . The grid 
structure used in our code is a multi-level set of properly 
nested, cell-centered uniform grids. We use a standard 
second-order finite difference stencil for the Laplacian op- 
erator and first order interpolation across the refinement 
level boundaries. The non-linearity of Eq. (p!9|) is ad- 
dressed by performing Newton- Raphson iterations. A 
brief description of our FMR implementation is summa- 
rized in Appendix [Aj to which we refer the interested 
reader. 



To check the consistency of solutions obtained with our 
FMR code we calculate the following diagnostic quanti- 
ties: 

The ADM mass is given by 

Madm = - ^ / 9'^^^^ j 

where we have applied Gauss' theorem to convert the 
surface integral into a volume integral. The actual ex- 
pression we use to calculate the ADM mass volume inte- 
gral is (|29|) with replaced by the right-hand-side of 

Eq. (uni). 

The Komar mass is given by 

Mk = ^ (f d'adSi = ^ [ V^aSx, (30) 
47r 47r J 

where again we have applied Gauss' theorem in the last 
step. By use of Eqs. (p!9|) and (|2T]) we find 

V^a = -^ViaV'^ + AiTa^f^{p + S). (31) 

The actual expression we use to calculate the Komar 
mass volume integral is (|3Q|) with V^a replaced by the 
right-hand-side of Eq. (f3T]) . 

Finally, the total baryon mass is given by ^] 



Mo = I poau^'^^d^x, 
Im 



(32) 



where A4 means that the integration is carried over the 
support of the matter. 



2. Code Testing 

Gauss' theorem constitutes a strong consistency check 
for our FMR code. To demonstrate that the solutions 
obtained with our elliptic code satisfy Gauss' theorem 
and achieve second-order convergence we performed the 
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FIG. 2: Mass - radius relationship of TOV stars generated 
by the piecewise polytropic EOS ([T5|). Plotted are curves cor- 
responding to three versions of the EOS where the ratio of 
the isotropic radius of a 0.98Mq pWD to that of a 1.5Mq 
NS is 20:1, 10:1 and 5:1. The corresponding EOS parame- 
ters are given in Tabled The open points correspond to the 
NS models studied in this paper, which have masses 1.4M0, 
1.5M0 and l.GM©. The solid points correspond to the pWD 
models considered in this work, which all have the same mass: 
O.98M0. 



A. Basic Equations 

The formulation and numerical scheme for our simula- 
tions are the same as those already reported in j22|, 
[67|, to which the reader may refer for details. Here we 
introduce our notation and summarize our method. 

We use the 3+1 formulation of general relativity where 
the metric is decomposed in the form of Eq. ([17]), and 
where the fundamental dynamical variables for the met- 
ric evolution are the spatial three-metric jij and extrin- 
sic curvature Kij. We adopt the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) formalism [11, 68, 69], in 
which the evolution variables are the conformal expo- 
nent (/) = ln(7)/12, the conformal 3-metric ^ij = e~^^jij, 
three auxiliary functions F* = — 7*-^ ,j, the trace of the 
extrinsic curvature and the trace- free part of the con- 



formal extrinsic curvature Aij 



7 = det(7ij). The full spacetime metric g^^ is re- 

gfiu^n^Uj,, where 



Here 

lated to the three-metric j^jy by j^jy 
the future-directed, timelike unit vector normal to the 
time slice can be written in terms of the lapse a and shift 
as = a~^(l, — The evolution equations of these 
BSSN variables are given by Eqs. (9)-(13) in [22]. 

We adopt standard puncture gauge conditions: an ad- 
vective "l+log" slicing condition for the lapse and a 
"Gamma- freezing" condition for the shift [70j. Thus, we 
have 



doa = —2aK , 
dop' = (3/4)5^ , 
doB' = doV-riB' 



(33) 
(34) 
(35) 



where do = dt — P^dj. We set the r] parameter to 
O.Olkm"^ for all simulations presented in this work. 

The fundamental matter variables are the rest-mass 
density po, specific internal energy e, pressure P, and 
four- velocity u^. We write the stress-energy tensor as 



following test. Employing the 10 : 1 piecewise poly- 
tropic EOS we constructed TOV NS solutions of vari- 
ous masses with a ID TOV integrator. We then used 
second-order polynomial interpolation to set up the rest- 
mass density profiles in our FMR elliptic code and solve 
Eqs. (fT9|) and ([2T]) for the conformal factor and lapse 
function. We set up grids with five levels of refine- 
ment {nl = 5) centered on the NS, and three differ- 
ent resolutions nx = ny = nz = 32, Ax^ = 2.2km, 
nx = ny = nz = 64, Ax^ = 1.1km and nx = ny = nz = 
128, Ax^ = 0.55km. In Fig. [3] we demonstrate that our 
numerical solutions satisfy Gauss' theorem, are in agree- 
ment with the TOV integration, and that our code is 
2nd-order convergent for this test. To generate the plot 
we used the results of the ADM mass integration but the 
convergence properties remain the same when we use the 
results of the Komar mass integration. 



(36) 



where h = 1 + e+P/po is the specific enthalpy and e is the 
total energy density. In our numerical implementation of 
the hydrodynamics equations, we evolve the "conserva- 
tive" variables p*, Si^ and f. They are defined as 



(37) 
(38) 
(39) 



The evolution equations for these variables are given by 
Eqs. (21)-(24) in 0. 

The EOS we adopt for the evolution has both a thermal 
and cold contribution, i.e.. 



P = Pth + Pc 



cold? 



(40) 
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FIG. 3: Convergence test for the FMR elliptic code using 
TOV stars. Five levels of refinement have been used for this 
test. Upper panel: Fractional difference between the volume 
integral Mvoi and surface integral Mgurf for the ADM mass, 
calculated with our FMR elliptic code, versus the gravita- 
tional mass Mtov calculated with a ID integrator of the TOV 
equations. The plot demonstrates satisfaction of Gauss' the- 
orem and second-order convergence. Lower panel: Fractional 
difference between Mvoi and Mtov versus Mtov- The plot 
demonstrates equality of Mvoi and Mtov and second-order 
convergence of our FMR elliptic code to the ID (exact) re- 
sult. In both panels three resolutions are plotted: low = 32"^, 
med = 64^, high = 128^. Resolutions 32^ and 64^ have been 
rescaled with a factor of 1/16 and 1/4 respectively, so that 
they overlap with resolution 128^. 



where Pcoid is given by Eq. ([15]) and the thermal pressure 
is given by 



Pth = (^th - l)po(e - ecoid), 



where 



Ccold 



I- 



cold 



d{l/po) 



(41) 



(42) 



We set Fth = 1-66 5/3) in all our simulations, i.e., we 
set it equal to the Fi exponent of the 10 : 1 EOS, appro- 
priate either for nonrelativistic cold degenerate electrons 
or (shock) heated, ideal nondegenerate baryons. Eq. (|4Q|) 
reduces to our piecewise polytropic law Eq. ([15]) for the 
initial (cold) NS and pWD matter. 



B. Evolution of the metric and hydrodynamics 

We evolve the BSSN equations using fourth-order ac- 
curate, centered finite-differencing stencils, except on 



shift advection terms, where fourth-order accurate up- 
wind stencils are applied. We apply Sommerfeld out- 
going wave boundary conditions on all BSSN fields, as 
in [22]. Our code is embedded in the Cactus paralleliza- 
tion framework [H], and our fourth-order Runge-Kutta 
timestepping is managed by the MoL (Method of Lines) 
thorn, with the CFL number set to 0.45 in ah pWDNS 
simulations. We use the Carpet [72| infrastructure to 
implement the moving-box adaptive mesh refinement. In 
all AMR simulations presented here, we use second-order 
temporal prolongation, coupled with fifth-order spatial 
prolongation, and impose equatorial symmetry to reduce 
the computational cost. 

We write the general relativistic hydrodynamics equa- 
tions in conservative form. They are evolved via a high- 
resolution shock-capturing (HRSC) technique 0, [63| 
that employs the piecewise parabolic (PPM) reconstruc- 
tion scheme [73*], coupled to the Harten, Lax, and van 
Leer (HLL) approximate Riemman solver [74]. The 
adopted hydrodynamic scheme is second-order accurate. 
To stabilize our hydrodynamic scheme in regions where 
there is no matter, a tenuous atmosphere is maintained 
on our grid, with a density floor patm set to 10~^^ 
times the initial maximum density on our grid. The ini- 
tial atmospheric pressure Patm is set by using the cold 
EOS ([T5]) . Throughout the evolution, we impose limits 
on the pressure to prevent spurious heating and nega- 
tive values of the internal energy e. Specifically, we re- 
quire Pmin < P < Pmax, whcrC P^ax = lOPcold and 
Pmin 

= O.SPcoid, where Pcoid is the pressure calculated 
using the cold EOS ([T5]) . Whenever P exceeds Pmax or 



drops below Pmin, we reset P to Pmax or Pmin, respec- 
tively. Following [75] we impose the upper pressure limits 
only in regions where the rest-mass density remains very 
low (po < lOOpatm), but wc imposc the lower limit every- 
where on our grid. 

At each timestep, the "primitive variables" po, P, and 
v'^ must be recovered from the "conservative" variables 
p^^ and Si. We perform the inversion numerically as 
specified in [58i]. We use the same technique as in [22| 
to ensure that the values of Si and f yield physically 
valid primitive variables, except we reset f to 10~^^fo,max 
(where fo,max is the maximum value of f initially) when 
either Si or f is unphysical [i.e., violate one of the inequal- 
ities (34) or (35) in [22[]. The restrictions are usually 
imposed only in the low-density atmosphere. 

It is instructive to discuss the mathematical structure 
of the system of hydrodynamic equations when a piece- 
wise polytropic EOS ([T5]) is used. According to |76|, [t^, 
when the fluxes of the conservation laws are non-smooth, 
split waves and composite structures may be present in 
the solutions. In these cases a numerical solution may 
not converge to the correct solution. In our case the 
fluxes are not smooth everywhere because EOS ([T5]) is 
non-smooth (it is continuous but not differntiable) at the 
turning points pi z = 1, 2. Away from the turning points 
the fluxes are smooth, therefore there may be some con- 
cern when non-linear waves cross the transition densities. 
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However, according to [77] our adopted numerical scheme 
should be able to handle such composite structures, if 
they ever arise, and our solutions should converge to the 
correct continuum solution. 

To study this effect we constructed an EOS which is 
similar to Eq. (p!5]) . but where the pressure discontinuities 
are "smoothed out" at the turning points pi^{i = 1,2), 
such that the EOS becomes a smooth, once (or twice) dif- 
ferentiable function of the rest-mass density. We perform 
such a smoothing operation using a cubic (or quintic) 
spline over a density interval [pi{l — e),p^(l + e)], where 
e > 0. We chose e to be sufficiently small so that the 
smoothed EOS mimics as closely as possible EOS (p!5]) . 
but large enough to avoid roundoff errors due to very 
large gradients. Setting up several generalized Riemman 
problems, we found that the numerical solutions obtained 
using EOS (dSj) converge to those obtained when using 
its smooth counterpart and the two can hardly be dis- 
tinguished for the resolutions considered. Therefore, our 
numerical schemes in conjunction with EOS ([15]) are able 
to capture the correct solution, in that they are almost 
identical to the solutions obtained with the smooth coun- 
terpart of (p!5|) . For the details of this analysis, we refer 
the interested reader to Appendix [Bl 



C. Evolution Diagnostics 

During the evolution, we monitor the normalized 
Hamiltonian and momentum constraints as defined in 
Eqs. (40)-(43) of [22^. 

We also monitor the ADM mass and angular mo- 
mentum of the system, which can be calculated dur- 
ing the evolution by surface integrals at a large distance 
(Eqs. (37) and (39) of [22]). The equations used to 
calculate the ADM mass and angular momentum with 
minimal numerical noise are as follows [llj 

-^X^Akmdntn + SwX^Sn] • (44) 

Here V is the volume within a distant surface, tjj = e^^ 
p = n^UyT^^ , Si = —Ujj^jijyT^^ , R is the Ricci scalar 
associated with jij , and Tijk are Christoffel symbols as- 
sociated with jij . 

In this work we only focus on head-on collisions, so 
there is no angular momentum involved. However, our 
simulations are three-dimensional, so there is no guaran- 
tee that Ji will remain 0. In order to quantify violations 



of Ji = we normalize the angular momentum, com- 
puted via dHI, with the angular momentum a pWDNS 
system would have, if the binary components were New- 
tonian point masses in circular orbit at the initial sepa- 
ration A, i.e., 

where the total mass Mt is taken to be the sum of ADM 
masses of the isolated stars. 

When hydrodynamic matter is evolved on a fixed uni- 
form grid, our hydrodynamic scheme guarantees that the 
rest mass Mq is conserved to machine roundoff error. 
This strict conservation is no longer maintained in an 
AMR grid, where spatial and temporal prolongation is 
performed at the refinement boundaries. Hence, we also 
monitor the rest mass 

Mo = p.d^x (46) 

during the evolution. Rest-mass conservation is also vio- 
lated whenever po is reset to the atmosphere value. This 
usually happens only in the very low-density atmosphere. 
The low-density regions do not affect rest-mass conser- 
vation significantly. 

Shocks occur when stars collide. We measure the 
entropy generated by shocks via the quantity K = 
P/Pcoid ^ 1, where Pcoid is the pressure associated with 
the cold EOS that characterizes the initial matter (see 
Eq. (dSD). 



VI. CASES AND RESULTS 

A. Initial configurations 

We perform a number of pWDNS head-on collision 
simulations varying the initial configurations, so that we 
can study the effect of the pWD compaction and NS mass 
on the final outcome separately. Table [III outlines the 
physical parameters for the cases considered in this work, 
and Table [IIll presents the AMR grid structure used in 
each case. 

To generate initial data for our cases we ffist choose the 
ADM mass Mns of the NS and the ADM mass Mwd of 
the pWD (see Table [III) solve the TOV equations for 
each star in isotropic coordinates to prepare the rest-mass 
density distribution for the NS and the pWD separately. 
We then use second-order polynomial interpolation to in- 
terpolate the rest-mass density profiles onto the nested 
grids of our FMR elliptic initial value code and solve 
Eqs. (p!9|) and ([2T]) for ^ and a. The two stars are placed 
at coordinate separation A and such that the Newtonian 
center of mass of the system is identified with the origin 
of the coordinate system. Once a solution is achieved by 
the FMR code for the initial metric, we map po, ^ and a 
from the elliptic code grids onto the evolution grids using 
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TABLE II: Summary of initial configurations. Mns (Mwd) stands for the ADM mass of an isolated NS (pWD)^"^ i?NS (i?WD) 
is the isotropic radius of an isolated NS (pWD), Cns (Cwd) is the compaction of an isolated NS (pWD), where the compaction 
is defined as the ratio of the ADM mass of the isolated star to its areal radius. Madm is the ADM mass of the system and 
A the initial binary separation in isotropic coordinates. All cases have exactly the same coordinate separation of 586.9km to 
allow for comparison. Cases Al, A2, A3 have been produced with the 10:1 EOS, while cases B and C have been produced with 
the 5:1 and 20:1 EOSs, respectively. 



Case 


Mns/M© 


Mwd /Mo 






Rwb/Rns 


i?WD/MADM 


Madm /Mo 


A/i?WD 


Al 


1.4 


0.98 


0.111 


0.010 


8.88 


41.18 


2.413 


4.000 


A2 


1.5 


0.98 


0.130 


0.010 


9.96 


39.36 


2.524 


4.000 


A3 


1.6 


0.98 


0.151 


0.010 


11.15 


37.46 


2.652 


4.000 


B 


1.5 


0.98 


0.130 


0.019 


4.99 


19.76 


2.524 


7.967 


C 


1.5 


0.98 


0.130 


0.005 


20.01 


79.08 


2.524 


1.991 



Here we list the ADM masses, isotropic radii and compactions of the isolated (TOV) stars whose rest-mass density profiles 
we used to generate initial data for ^ and a for a given case. 



TABLE III: Grid configurations used in our simulations. Here M is the sum of the ADM masses of the isolated stars. Max. res. 
is the grid spacing in the innermost refinement box surrounding the NS, A^ns denotes the number of grid points covering the 
diameter of the NS initially, and A^vD denotes the number of grid points covering the diameter of the pWD initially. The 
smallest outer boundary distance corresponds to case A3 and is 1028M. 



Case 


M/Mo 


Grid Hierarchy (in units of M)^^) 


Max. res. 




Al 


2.38 


(534.33, 267.16, 133.58, 66.79, 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) 


M/6.71 


63 


35 


A2a 


2.48 


(510.62, 255.31, 127.65 , 63.83, 34.81[N/A], 18.86[N/A], 10.15[N/A], 6.890[N/A]) 


M/5.52 


44 


28 


A2b 


2.48 


(534.33, 267.16, 133.58, 66.79, 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) 


M/6.71 


56 


35 


A3 


2.58 


(467.27, 233.64,116.82, 58.41, 29.20 [N/A], 15.58[N/A], 8.518[N/A], 5.841[N/A]) 


M/8.22 


56 


38 


B 


2.48 


(534.33, 267.16, 133.58, 66.79, 35.78[31.93], 19.08[N/A], 10.44[N/A], 7.156[N/A]) 


M/6.71 


56 


35 


C 


2.48 


(534.33, 267.16, 133.58, 66.79[N/A], 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) 


M/6.71 


56 


35 



There are two sets of nested refinement boxes: one centered on the NS and one on the pWD. This column specifies the 
half side length of the refinement boxes centered on both the NS and pWD. When the side length around the pWD is 
different, we specify the pWD half side length in square brackets. When there is no corresponding pWD refinement box (as 
the pWD is much larger than the NS), we write [N/A] for that box. 



second-order polynomial interpolation. We always make 
sure that the resolution of the initial data grids is higher 
than the resolution of the evolution grids. The surfaces of 
the stars are a locus of rapidly decreasing density gradi- 
ents. As a result, small oscillations due to interpolation 
may arise and lead to negative rest-mass density. To 
cure this, we set the density po equal to the tenuous at- 
mosphere density that we maintain on the grid whenever 
\po/pi,c\ < 10~^^, where po is the value of the density af- 
ter the interpolation and p^^c, ^ = WD, NS, is the central 
density of the WD or NS. We do not find such oscillations 
when interpolating the gravitational fields, which is most 
likely due to the fact that these are sufficiently smooth. 

In all cases we require that the sum of the ADM masses 
of the isolated stars be larger than the maximum gravi- 
tational mass of I.8M0 that our cold EOS can support 
(see Fig.[2j). There exist at least 18 observed WDNS sys- 
tems that satisfy this requirement. Since typical NS and 
WD masses in massive WDNS binaries lie in the range 
I.3M0 - I.6M0 and O.5M0 - I.IM0 respectively (see 
Sec. H]), we choose pWD rest-mass density profiles that 



correspond to an ADM mass of O.98M0 in isolation and 
keep it fixed in all cases we study. This almost fixes 
the pWD rest mass, because of the small compaction 
(< 0.02) of the pWDs we consider. The pWD rest-mass 
variation from case B to case C, due to fixing the ADM 
mass, is 0.7%. We vary only the pWD compaction, i.e., 
the EOS, and the NS mass. The reason why we chose the 
pWD mass to be O.98M0 is that the ratio of the isotropic 
radius of a pWD of such mass to the isotropic radius of 
a I.5M0 NS is ^ 10 for the 10:1 EOS. 

Another quantity that we fix in all our simulations is 
the initial coordinate separation A of the two compo- 
nents. This almost fixes the kinetic energy of the stars 
when they collide. In particular, we set A = 4i?wD,A- 
Here Rwb^a denotes the isotropic radius of the spheri- 
cal O.98M0 pWD used in cases Al, A2, A3. We choose 
the initial separation this way because we want the stars 
to be sufficiently far apart so that spherical TOV ini- 
tial models remain in near equilibrium, and at the same 
time, simulate the collision within reasonable time scales, 
as the head-on collision time scale varies as ~ 

^3/2 ( 

see 
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FIG. 4: Snapshots of rest-mass density profiles at selected times for case A2. The contours represent the rest-mass density 
in the orbital plane, plotted according to po — po,maxlO~°"^^-^~°"^^ {j = 0, 1, ... ,9). The color sequence dark red, red, orange, 
yellow, green, light green, blue and light blue implies a sequence from higher to lower values. This roughly corresponds to 
darker greyscaling for higher values. The maximum initial NS density is /?o,max = 4.6454/>nuc- The last two snapshots are near 
the end of the simulation, and they demonstrate that the density contours within a radius of about 150km remain unchanged. 
Here M = 2ASMq = 3.662km = 1.222 x 10~^s is the sum of the ADM masses of the isolated stars. 



Eq. ©). 

li A = 4i?wD,A, the NS tidal field in the vicinity of 
the pWD is small, validating our assumption that an 
equilibrium pWD is nearly spherical. To see this, let us 
calculate the ratio of the tidal force of the NS on the 
surface of the pWD to the surface gravity of the pWD 



5%, 



(47) 



where we used Mns = l.SM©, Mwd = l^o, A = 
3i?WD,A- Hence, any deviations from sphericity should 
be small. The assumption of sphericity for case B is 
even better because in this case A 8-Rvvd b ^ but worse 
for case C, where A ^ 2i^wD,c- In principle, we could 
increase the separation so that the sphericity approxima- 
tion becomes better for all cases, but if the final remnant 
does not collapse promptly to form a BH starting at close 
separations, it is unlikely that it will collapse if the initial 
separation is larger. This is due to the fact that for larger 
separations the kinetic energy at collision will be larger. 



generating more shock heating that will work to prevent 
prompt collapse. 

To summarize, the set of cases Al, A2 and A3 probe 
the effect of the NS mass on the final outcome, whereas 
the set of cases A2, B and C probe the effect of the pWD 
compaction on the final outcome. In the following sec- 
tions we summarize the results of our simulations. 



B. Dynamics of collision and effects of the NS mass 

Here we describe the effects of the NS mass on the 
dynamics of pWDNS head-on collisions. We find that 
about 18% of the initial total mass escapes to infinity for 
all cases Al, A2 and A3. Nevertheless, the initial total 
mass in these cases is large enough to guarantee that the 
final total mass of the pWDNS remnant still exceeds the 
maximum mass that our cold EOS can support. How- 
ever, prompt collapse to a black hole does not take place 
in any of the cases studied because strong shock heat- 
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FIG. 5: Normalized angular momentum vs time. Jz and Jz,c 
are given by Eqs. (|44)) and (|^ . respectively. The solid, 
dashed and dotted lines correspond to cases A2, B and C, 
respectively. Here M = 2A8Mq = 3.662km = 1.222 x 10"^s. 



ing gives rise to a hot remnant. Ultimate collapse to a 
BH is almost certain after the remnant has cooled. The 
outcome of the three cases Al, A2, A3 is a TZIO. 

Overall, cases Al, A2 and A3 are qualitatively similar 
and for this reason we mainly describe case A2 as repre- 
sentative of this class of our simulations. Furthermore, 
our study of case A2 with two different resolutions (see 
Table HITl) shows the results to be qualitatively insensitive 
to resolution indicating that the resolutions used in our 
simulations are sufficiently high. In what follows all case 
A2 plots correspond to the high resolution run of case 
A2, i.e., case A2b. 

In general, the head-on collision of pWDNS systems 
can be decomposed into three phases: the acceleration, 
the plunge, and the quasiequilibrium phase. 

During the acceleration phase, the two stars accelerate 
toward one another starting from rest. The separation 
decreases as a function of time and this phase ends when 
the two stars first make contact. 

As the separation decreases, the increasing NS tidal 
field strongly distorts the pWD. This can be seen in the 
equatorial rest-mass density contours of Fig. |4l In the 
insets of Fig. 21 the NS interior is almost unchanged dur- 
ing this phase. In reality, it oscillates but is not tidally 
distorted by the pWD. Nevertheless, the NS atmosphere 
does expand. The insets also show that due to numerical 
errors the NS veers slightly off the x-axis, which is the 
collision axis in our simulations. In general, in all our 
simulations both the NS and the pWD wiggle around 



the X-axis. The amplitude of the NS wiggling motion is 
at most 1% of the pWD radius, while the amplitude of 
the pWD off-axis motion is less than 0.1% of its radius. 
Hence, the collision is practically head-on. 

It is likely that this lack of symmetry is due to small 
asymmetries introduced when mapping the initial data 
onto the evolution grids via 2nd order interpolation. This 
is because 2nd order interpolation requires the use of 3 
grid points (per direction) that surround the point to 
which one interpolates. However, the effect is small and 
our results cannot change qualitatively due to this small 
asymmetry. 

Along with this small off-axis motion, the pWDNS sys- 
tem acquires a small amount of spurious angular momen- 
tum. Fig. [5] shows the normalized z component of the 
angular momentum of the system and demonstrates that 
it is always less than 3% (see Eq. (|45]) ). In addition to 
conserving the angular momentum to within 3%, the nor- 
malized Hamiltonian constraint violations remain smaller 
than 1% and the normalized momentum constraint vio- 
lations smaller than 3%. These results hold for all cases 
studied in this work. 

During the plunge phase the NS penetrates the pWD, 
launching strong shocks that sweep through and heat 
the interior of the pWD. The NS outermost layers are 
stripped when they encounter the dense central parts of 
the pWD, and the NS is compressed. We find that at 
maximum compression in case A2, the NS central density 
only increases by about 8% of the initial central density. 

Eventually, strong shocks sweep through the entire 
pWD interior and then transfer linear momentum to the 
pWD outer layers, a large fraction of which receives suf- 
ficient momentum to escape to infinity. This can be seen 
in Fig. [6l where a snapshot is shown of the total energy 
per unit mass U = —uq — I (subtracting the rest-mass en- 
ergy) on the equator long after the collision, when ejected 
material has already started crossing the outer boundary 
of the computational domain. Unbound {U > 0) mat- 
ter covers most of the computational domain, as shown 
in Fig. [6l The rest-mass density of the ejected material 
is of order lO^^pnuc, but the total mass that escapes to 
infinity is large. This is demonstrated in Fig. [71 which 
shows the fractional change in the rest mass as a function 
of time. We find the amount of matter that escapes in 
cases Al, A2 and A3 is c:^ 18% of the initial rest mass 
when we extrapolate our results to late times. 

The thermal energy deposited into the ejected material 
is significant, with K = P/Pcoid > 40. As the ejected 
matter comes from the WD outermost layers, its den- 
sity is very low. This implies that its initial pre-shocked 
sound speed is small. As a result, the Mach number of 
the ejected material can be very large prior to shock heat- 
ing, and so shock heating is very strong, i.e., K increases 
from 1 initially to greater than 40 (see also discussion in 
Appendix B of [23]). 

The time scale for shock heating to equilibrate must 
be of order few times the dynamical (free-fall) time scale 
of the WD (see Eq. ([2])), as this is the only relevant 
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FIG. 6: Snapshot of total energy per unit mass U = —uq — 1 (subtracting the rest-mass energy) for cases A2 and B, when 
matter has already started crossing the outer boundary of the computational domain. Matter with [/ < is bound, while 
U > implies that matter is unbound. The contours represent the total energy per unit mass in the equatorial plane, plotted 
according to [/ = UminlO^'^^^ (j = 0, 1, . . . , 9). The color code here is the same as that defined in Fig. [H The white spaces in 
the centers of the plots correspond to bound matter (U < 0). We chose the cutoff value Umin = 10~^. The size of the bound 
matter area is insensitive to the choice of Umin- Here M — 2.48M0 = 3.662km = 1.222 x 10~^s. 



timescale. For cases A, the WD dynamical timescale 
is roughly 400M. Our computations show that it actu- 
ally takes about 800M-1000M for the star to equilibrate, 
which is consistent with the estimate above. 

Material that did not receive sufficient thrust to escape 
to infinity starts to rain down onto the NS and pWD 
matter and the plunge phase ends when this process is 
over. 

During the quasiequilibrium phase the remnant settles 
into a spherical quasiequilibrium object whose outer lay- 
ers oscillate. This can be seen in the two final snapshots 
of Fig. m where we show that the inner equatorial rest- 
mass density contours do not change with time, while 
the outer layers change only a little. Fig. [8] plots xy, xz 
and yz rest-mass density contours. Notice that in the 
adopted gauge, the remnant appears to be spherical. 

The pWDNS final remnant consists of a cold NS core 
with a hot mantle on top. This is demonstrated by 
the plots in the last row of Fig. [H where we quan- 
tify the results of shock heating by showing contours of 
K = P/Pcoid- Within a radius of 100 km from the center 
of mass of the remnant, K ranges from unity to about 
15 for case A2. In all cases 1 at the center of the 

remnant, while it becomes larger than unity as we move 
outwards from the center. We refer to this spherical con- 
figuration as a Thorne-Zytkow-like object. 

Even though a large fraction of the initial mass escapes, 
the final total rest mass well exceeds the maximum rest 
mass of 1.92M0 our cold EOS can support. In Fig. [9] 



we show the rest mass of the remnant as a function of 
time and for various spatial domains. This figure demon- 
strates that the rest mass within a radius of 220km ac- 
counts for more than 90% of the final total rest mass and 
is greater than 1.92M0. However, the pWDNS remnant 
does not collapse promptly to form a black hole, because 
of extra support provided thermal pressure. Delayed col- 
lapse to a black hole is almost certain after the pWDNS 
remnant has cooled. 

Finally, we note that it has been suggested in [tS*, |79| 
that GWs may arise from shocks. Even though the dis- 
cussion in these studies focused on core collapse super- 
novae, the appearance of strong shocks in our case can 
also generate GWs. However, here we do not calculate 
the GW signature because what is really interesting from 
an observational and astrophysical point of view is the 
GW signature in the circular binary WDNS case, not 
the head-on case we consider. In [33j we did calculate 
GWs from the inspiral phase of binary WDNS systems. 
General relativistic computations of the merger of circu- 
lar binary WDNSs will be the subject of a forthcoming 
paper. 



C. Effects of the pWD compaction 

Here we describe the effects of the pWD compaction 
on the dynamics of pWDNS head-on collisions. Overall, 
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TABLE IV: Summary of pWD compaction study. Here Cwd is the compaction of an isolated pWD (see Table . K — P/Pcoid 
at the end of simulations^"-*, Tp is the peak temperature at the end of simulations, Mo(0) is the initial total rest mass, 
AMo = |Mo,f — Mo (0)1, where Mo,f is the final total rest mass, Mo,r<22o is the mass enclosed within 220 km from the center of 
mass of the remnant at the end of the simulations, po,c is the final value of the central rest-mass density *^ , and ttmin the final 
value of the minimum of the lapse function. 



Case 




K 


Tp(10^^ °K)(^) 


AMo /Mo (0) 


Mo,r<22o/M0 


PO,c/ Pnuc 


Q^min 


B 


0.0190 


[1,35] 


3.7 


14% 


2.180 


4.91 


0.570 


A2 


0.0100 


[1,15] 


3.2 


18% 


2.035 


4.49 


0.595 


C 


0.0049 


[1,10] 


3.0 


18% 


1.900 


4.10 


0.609 



The K column lists the range of values which K obtains within a radius of 100 km from the centers of mass of the 
remnants. 

Pnuc = 2 X IQl^g/cm^. 

(c) por realistic WDNS collisions we expect Tp ~ 10® °K (see discussion following Eq. (ET))). 
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FIG. 7: Fraction of rest mass lost vs time. Here AMo = 
I Mo — Mo (0)1, where Mo(0) is the initial total rest mass. 
Small changes in the rest mass until approximately 5000M 
for cases Al, A2, A3, B and lOOOOM for case C are due to 
interpolations when matter crosses refinement levels and in- 
accuracies in evolving the very low-density atmosphere. At 
the end of the simulations the amount of mass ejected is 
13.4% of the initial rest mass in case B and ranges from 
16.1%- 16.7% of the initial rest mass for the other cases. 
Extrapolating the results to late t we find that in case B 
AMo /Mo (0) ^ 14%, AMo /Mo (0) 18% in ah other cases. 



Here M = 2.48M( 







3.662km = 1.222 x lO-^'s. 



our findings are qualitatively similar to those of case A2 
described in the previous section. An appreciable frac- 
tion of the initial total mass escapes to infinity, but the 
final total mass of the pWDNS remnant still exceeds the 
maximum mass that our cold EOS can support. Prompt 
collapse to a black hole does not take place either in case 



B or in case C, because strong shock heating gives rise to 
a hot remnant. The outcome of cases B and C is again a 
TZIO. 

The three phases of the head-on collision we described 
in the previous section apply here, too. For this reason 
we now focus our discussion on describing the differences 
between cases B, A2 and C, i.e., in order of decreasing 
pWD compaction. 

The tidal acceleration, which the pWD experiences due 
to the NS tidal field, increases as the pWD compaction 
decreases. This is because the initial coordinate separa- 
tion is the same for cases B, A2 and C. As a result, the 
acceleration phase is shorter for larger pWD compaction. 
This can be seen in Figs. [10] and [TTJ where equatorial 
rest-mass density contours are plotted. 

Shock heating far from the core of the remnants, as 
quantified by is somewhat less intense as the pWD 
compaction decreases. The shorter acceleration phase 
implies that the relative speed of the two components 
at the beginning of the plunge phase is a little smaller, 
which in turn leads to weaker shocks. 

During the plunge phase, the NS interior is less affected 
by decreasing pWD compaction. This can be seen (a) in 
the insets in Figs. [10] and [TTJ where in case C the post- 
plunge structure of the NS core is almost the same as that 
showed in the pre-plunge snapshots, while this is not true 
for case B, and (b) by the variation in the NS central 
density (po,c); in particular, we find that at maximum 
compression the NS central density increases by about 
42% in case B, 8% in case A2 and 5% in case C. These 
results can be easily interpreted because in a sequence 
of pWDNS head-on collisions where the NS is fixed and 
the size of the pWD increases with fixed mass, the NS 
gradually encounters thinner and thinner material, and 
hence changes to the NS structure become less and less 
important. 

Were the system mass loss to vary appreciably with 
pWD size, we might expect a corresponding variation in 
po,c- However, such a mass loss variation is not observed, 
as we discuss next. 

As in case A2, in both cases B and C, a large frac- 
tion of the initial mass eventually escapes to infinity (see 
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FIG. 8: First three rows: Snapshots of rest-mass density profiles at the end of the simulation for cases A2, B, C. The contours 
represent the rest-mass density in the YZ plane (first row), in the XZ plane (second row) and in the XY plane (third row) 
plotted according to po = Po,maxl0~'^"^^'^~'^"^^ (j = 0, 1, . . . , 9). The maximum initial NS density is po,max = 4.6454/>nuc- These 
snapshots demonstrate that in the adopted gauge, the final object is roughly spherical. Last row: Snapshots of K = P/Pcoid 
profiles at the end of the simulation for cases A2, B, C. The contours represent K in the XY plane plotted according to 
K = IQO-^^^-^ (j = 0, 1, . . . , 9). It is evident that the core of the remnant remains cold {K c^i 1). K becomes larger than unity 
as we move outwards from the center of the objects, and shock heating is more intense in case B and less intense in case C. 
The color code used is the same as that defined in Fig. H Here M = 2A8Mq = 3.662 km = 1.222 x 10~^ s. 
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FIG. 9: Post-merger rest mass as a function of time. Here Mo, tot is the total rest mass in the entire computational domain and 
Mo,r<ro stands for the rest mass contained within a coordinate sphere of radius ro in units of km, centered on the remnant's 
center of mass. In all cases Mo,r<22o accounts for more than 90% of the final total rest mass and it is always greater than 
1.92M0 - the maximum rest mass that our cold EOS can support. Here M = 2.4&Mq = 3.662km = 1.222 x 10"^s. 



Fig. [6] for case B), but we do not find strong variations 
in the mass lost among cases B, A2 and C. We find that 
the amount of matter that escapes in case B is 14% of 
the initial rest mass, while case C loses 18%, (case A2 
loses 18%) (see Fig. [7]) of the initial rest mass. Given 
our earlier discussion that shocks in case B are stronger 
than those in case A2, and in turn shocks in case A2 are 
stronger than those in case C, this last result may sound 
contradictory, because one might expect that stronger 
shocks would eject more matter to infinity. The appar- 
ent contradiction can be resolved, if one considers that 
as the pWD compaction decreases, the pWD outer layers 
become less and less bound, and hence, it requires less 
energy to eject them to infinity. 

As in case A2, the remnants in cases B and C even- 
tually settle into spherical quasiequilibrium objects with 
oscillating outer layers (see Figs. fTOl [TT ]) . The sphericity 
of the remnants (in the adopted gauge) in cases B and C 
is demonstrated by the xy, xz and yz rest-mass density 
contours shown in Fig. [51 

The pWDNS remnants in both case B and case C con- 
sist of a cold NS core with a hot mantle on top. Thus, 
all cases lead to the formation of a TZIO. This is again 
demonstrated in the last row of Fig. [51 where contours of 
K = P/Pcoid are shown. Within a radius of about 100 km 
from the center of mass of the final remnants, G [1, 35] 
in case B, and K e [1, 10] in case C {K e [1, 15] in case 
A2). Thus, shock heating is overall strongest in case B, 
weaker in case A2 and even weaker in case C. 

In Fig. [9l the remnant rest masses in cases B and C are 
plotted as functions of time. The figure demonstrates 
that in both cases, the rest masses within a radius of 
220km account for more than 90% of the final total rest 
masses in both case B and case C, and are greater than 
1.92M0, i.e., the maximum rest mass that our cold EOS 
can support. The pWDNS remnant does not collapse 
promptly to form a black hole, because of the extra ther- 



mal pressure support. However, delayed collapse to a 
black hole is almost certain after the pWDNS remnant 
has cooled. 

Another important feature of Fig. [9l is that the amount 
of mass contained within a given radius from the center 
of mass of the remnant is larger for smaller initial pWDs. 
For example, within a radius of 220km the remnant mass 
is 2.I8M0 in case B, 2.O35M0 in case A2, and I.9OM0 
in case C. This in turn indicates that the higher the ini- 
tial pWD compaction the higher the core densities of the 
pWDNS remnant. This is supported by the rest-mass 
density contours shown in Fig. [51 and by the values of the 
final central rest-mass density. In particular, we find that 
the final central rest-mass density is 4.10pnuc in case C, 
4.49pnuc in case A2, and 4.91pnuc in case B. Thus, there 
is a variation in the final central density of 9.2% from 
case B to case A2, and 9.5% from case A2 to case C. 
Finally, it is also worth noting that the final minimum 
value of the lapse function, which is a good indicator of 
collapse, increases with increasing initial pWD size, too. 
Specifically, we find this value to be 0.57 in case B, 0.595 
in case A2, and 0.609 in case C. All these facts seem to 
indicate that as the pWD size increases towards realistic 
WD sizes the less likely it is for the pWDNS remnant to 
collapse. To demonstrate this trend more clearly we com- 
pile all aforementioned physical parameters of the final 
configurations in cases B, A2, and C in Table [TVl 

Hence, given the consistency in the results of cases 
B, A2 and C, i.e., the sequence of increasing pWD size 
with fixed pWD mass, we expect that as the parameters 
of our EOS vary, so as to describe realistic WDs, the 
result of the head-on collision of a massive WDNS system 
will most likely lead to formation of a quasiequilibrium 
TZIO. If the initial total mass of the system exceeds the 
maximum mass a cold EOS can support, then the TZIO 
may have mass that exceeds the maximum mass a cold 
EOS can support, so it will eventually undergo collapse 
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FIG. 10: Snapshots of rest-mass density contours at selected times for case B. The contours are plotted in the orbital plane, 
according to po = po,maxlO~°'^^-^~°'^^ (j = 0, 1, . . . , 9). The color code used is the same as that defined in Fig.U] The maximum 
initial NS density is po,max = 4.6454y9nuc- The last two snapshots, which take place near the end of the simulation, demonstrate 
that the density contours within a radius of about 150km remain unchanged. Here M = 2A8Mq = 3.662km = 1.222 x 10~^s 
is the sum of the isolated stars' ADM masses. 



to a black hole, but only after the remnant has cooled. 

To identify the dominant cooling mechanisms and/or 
relevant nuclear reaction networks, one would need to 
estimate the temperatures of these TZlOs. We can do 
this as follows. 

Using Eq. (|^Q|) and the definition of K we can calculate 
the specific thermal energy as 



{K - l)Pcold 

{Tth - l)po " 



(48) 



To estimate the temperature T of matter, we assume 
that we can model the temperature dependence of eth as 



3kT 
2mn 



J ; 

PO 



(49) 



where rrin is the mass of a nucleon, k is the Boltz- 
mann constant and a is the radiation constant. The 
first term represents the approximate thermal energy 
of the nucleons, and the second term accounts for the 
thermal energy due to relativist ic particles. The fac- 
tor / reflects the number of species of ultrarelativistic 



particles that contribute to the thermal energy. When 
T « 2me//c ~ lO^^K, where ttIq is the electron mass, 
thermal radiation is dominated by photons and / = 1. 
When T » 2me//c, electrons and positrons become 
ultrarelativistic and also contribute to radiation, and 
/ = 1 + 2 X (7/8) = 11/4. At sufficiently high tem- 
peratures and densities (T ^ 10^^K,po ^ lO^^g cm~^), 
neutrinos are generated copiously and become trapped, 
so, taking into account three flavors of neutrinos and an- 
tineutrinos, / = 11/4 + 6 x (7/8) = 8. 

In Fig. [12] we show the temperature profiles of the rem- 
nants in cases B, A2 and C, where it is clear that typical 
temperatures of our TZlOs are of order 10^^ ^K. This is 
expected as the total energy available for shock heating 
should be of order Mns^wd/^wd, i-e., the gravitational 
interaction energy when the two stars first make contact. 
The total thermal energy, Eth^ is then 

(Mns + Mwd),^ MnsMwd 
Eth ^ kT — . (50) 

From this last equation we can estimate the characteristic 
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FIG. 11: Rest-mass density contours in the orbital plane at selected times for case C. The contours are plotted according to 
po — po,maxlO~°"^^-^~°"^^ (j = 0, 1, . . . , 9). The color code used is the same as that defined in Fig. U] The maximum initial NS 
density is y9o,max = 4.6454/>nuc- The last two snapshots, which take place near the end of the simulation, demonstrate that the 
density contours within a radius of about 150km remain unchanged. Here M — 2.48M0 = 3.662km = 1.222 x 10~^s is the sum 
of the isolated stars' ADM masses. 



temperature as 



VII. SUMMARY AND DISCUSSION 



(51) 



Thus, all things being equal (no mass loss, same kinetic 
energy at collision, etc.) characteristic TZIO tempera- 
tures should be directly proportional to the pWD com- 
paction. Taking case A2 as an example, C 0.01 and 
q 2/3^ we find T 6.5 x 10^^ ^K, in rough agreement 
with our simulations. 

Using this scaling argument we can extrapolate our 
results to realistic WDNS head-on collisions. For a 
solar-mass WD obeying the Chandrasekhar EOS Cwd — 
3 X 10~^. Hence, we predict that typical temperatures 
in realistic head-on collisions of massive WDNS systems 
would be of order 10^ ^K. 



In this work we studied the dynamics of the head-on 
collision of WDNS binaries in full general relativity, aided 
by simulations that employ the Illinois AMR relativistic 
hydrodynamics code [23^, 58] . This study serves as a pre- 
lude to the circular binary WDNS problem which will be 
the subject of a forthcoming paper. 

Our primary focus is on compact objects whose total 
mass exceeds the maximum mass that a cold EOS can 
support and our goal is to determine whether a WDNS 
collision leads either to prompt collapse to a black hole 
or the formation of a quasi-equilibrium configuration that 
resembles a Thorne-Zytkow object (TZO) [57|, which we 
call Thorne-Zytkow-like object (TZIO). By a TZIO we 
mean a NS sitting at the center of a hot gaseous mantle 
composed of WD debris. 

Due to the vast range of dynamical time scales and 
length scales involved in this problem, realistic WDNS 
simulations (head-on or otherwise) are computationally 
prohibitive, if one employs current numerical relativity 
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FIG. 12: Temperature (T) profiles for cases B, A2 and C. The 
temperature is in units of 10^^ "^K. The profiles correspond 
to the values of T at the end of the simulations and along the 
X-axis (for x > Xc, where Xc is the x position of the center of 
mass of the remnant in each case) . It is clear that typical tem- 
peratures are of order 10^^ '^K. For realistic WDNS collisions 
we expect T ^ 10^ °K (see discussion following Eq. (|5ip ). 

techniques and available computational resources. For 
this reason, we tackle the problem using a different ap- 
proach. 

In particular, we constructed a piecewise polytropic 
EOS which captures the main physical features of NSs 
and, at the same time, scales down the size of a WD. 
We call these scaled-down models of WDs pseudo-WDs 
(pWDs). Using these pWDs, we can reduce the range 
of length scales and time scales involved, rendering the 
computations feasible. 

A pWD is not a realistic model of a WD. However, 
with our proposed parametrized EOS we can construct a 
sequence of pWD models with gradually increasing size 
and perform simulations that approach the realistic case. 
Then we can make predictions about the realistic case 
by extrapolating the results from this sequence of simu- 
lations. 

Using pWDs, we performed two sets of simulations. 
One set of our simulations studied the effects of the NS 
mass on the final outcome, when the pWD is kept fixed at 
a mass of O.QSM© and its size fixed at 146km. We choose 
three masses for the NS, namely 1.4M0, l.SM©, 1.6Mq 
(cases Al, A2 and A3, respectively). The other set of 
simulations studied the effect of the pWD compaction on 
the final outcome, when the NS is kept fixed at a mass of 
1.5M0. In the latter set of calculations, we choose three 
values for the ratio of the pWD to the NS radius, namely 



5:1, 10 : 1 and 20 : 1 (cases B, A2 and C, respectively). 

In general, the head-on collision of pWDNS systems 
can be decomposed in three phases: i) acceleration, ii) 
plunge and iii) quasiequilibrium. 

During the acceleration phase the two stars accelerate 
towards one another starting from rest. As the separation 
decreases the NS tidal field becomes so strong that the 
pWD becomes highly distorted, while the NS interior is 
almost unchanged. This phase ends when the NS and 
pWD first make contact. 

During the plunge phase the NS penetrates the pWD, 
launching strong shocks that sweep through and heat 
the interior of the pWD. The NS outermost layers are 
stripped after encountering the dense central parts of the 
pWD (see Fig. |4]), but the NS core is mostly unaffected, 
except when the compaction of the pWD is high (see 
Sec. IVI C[) . We find that the strong shocks sweeping 
the pWD transfer linear momentum to the outer pWD 
layers, causing a large amount of pWD matter to escape 
to infinity. In all calculations, we find the rest mass loss 
to be between 14% — 18% of the initial total rest mass. 
Material that did not escape to infinity accretes onto the 
underlying NS and pWD matter. 

Finally, during the quasiequilibrium phase, the rem- 
nant settles into a spherical quasiequilibrium object 
whose outermost layers undergo damped oscillations. 

Although a large fraction of the initial mass escapes, 
the final total rest mass still exceeds the maximum rest 
mass of 1.92M0 that our cold EOS can support (see 
Fig. [9]). However, the pWDNS remnant cannot collapse 
promptly to form a black hole, because it is hot. This re- 
sult is the same in all our simulations. However, we point 
out that delayed collapse to a black hole is almost cer- 
tain after the pWDNS remnant has cooled, but this will 
occur on a timescale much larger than a hydrodynamical 
timescale. 

The final object consists of a cold NS core surrounded 
by a hot mantle. We quantified the results of shock heat- 
ing by the ratio of the total pressure to the cold pressure 
K = P/Pcoid- In all cases K 1 at the center of the 
remnant, and becomes larger than unity away from the 
center. We refer to this nearly-spherical configuration as 
a Thorne-Zytkow-like object, and find this object at the 
end of all simulations, regardless of NS mass and pWD 
compaction. We find that within a radius of 100 km from 
the centers of mass of the remnants, K lies in the range 
[1, 15] in cases Al, A2 and A3, [1, 35] in case B and [1, 10] 
in case C. Using a simple model for the temperature de- 
pendence of the specific thermal energy we estimate the 
characteristic temperature of these objects to be of order 
10^^ ^K. Using a simple scaling argument (see Eq. (|5T]) ) 
we find that TZIO temperatures should be proportional 
to the compaction of the original pWD, so that in realistic 
WDNS head-on collisions typical remnant temperatures 
would be of order 10^ ^K. 

Furthermore, we find that the smaller the initial pWD 
compaction the smaller the core densities of the pWDNS 
remnant. This is supported by the rest-mass density con- 
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tours shown in Fig.[8]and by the values of the final central 
rest-mass density. In particular, we find that the final 
central rest-mass density decreases by 9.2% from case B 
to case A2, and 9.5% from case A2 to case C. In addition, 
the final minimum value of the lapse function, which is a 
good indicator of collapse, increases with increasing ini- 
tial pWD size, too. Specifically, we find this value to be 
0.57 in case B, 0.595 in case A2, and 0.609 in case C (see 
Table IIVI for a summary of physical parameters of the 
final configurations in cases B, A2, C). All these facts 
seem to indicate that as the pWD size increases towards 
realistic WD sizes the less likely it is for the pWDNS 
remnant to collapse. 

An important concern regards the invariance of these 
results with respect to larger initial separations. To fully 
resolve this, one would need to extend the simulations to 
wider separations, but this extension is outside the scope 
of the current work. However, this work gives us some 
qualitative idea about what might happen with larger ini- 
tial separations. Larger separations imply larger kinetic 
energies during the plunge phase, which in turn imply 
stronger shocks. Stronger shocks likely lead to larger 
mass loss and more intense shock heating. Therefore, our 
expectation is that head-on collisions of pWDNS systems 
starting at larger separations will also result in the for- 
mation of TZlOs and that such collisions would not lead 
to prompt formation of a black hole. 

Given the consistency in the results of cases B, A2 
and C, we expect that as the parameters of our EOS are 
adjusted such that pWDs more closely resemble realis- 
tic WDs, WDNS head-on collisions are likely to form a 
quasiequilibrium TZIO. If the initial total mass of the 
system well exceeds the maximum mass that a cold EOS 
can support, then the TZIO will most likely have mass ex- 
ceeding the maximum mass supportable by a cold EOS, 
eventually collapsing to a black hole after the remnant 
has cooled. 

We conclude by stressing that we cannot use the re- 
sults of this work to make definite predictions about ei- 
ther the pWDNS or the realistic WDNS circular binary 
problem. One might speculate that shock heating will 
be minimized in such a scenario, and hence it may result 
in prompt collapse to a black hole. However, sufficient 
angular momentum must be shed in the circular binary 
case in order for the object to promptly form a black 
hole. To resolve these issues, hydrodynamic simulations 
in full general relativity must be performed and will be 
the focus of a forthcoming paper. 
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Appendix A: Initial data Code description 

In this appendix we describe some details of the fixed 
mesh refinement (FMR), finite difference code we devel- 
oped for generating general relativistic initial data. 

The grid structure we use is a multi-level set of prop- 
erly nested, cell-centered uniform grids. Each grid cor- 
responds to one level of refinement labeled by the level 
number il = (0, 1, 2, . . . , n/ — 1), where nl is the total 
number of levels, il = corresponds to the coarsest level 
and il = nl — 1 to the finest one. All levels have the same 
number of grid points nx, ny^ G Z in the x, y and z 
directions respectively. The coordinates on our grid are 
defined as follows 



z = 0, 1, . . . , nx — 1, 
j = 0, 1, . . .,n?/ - 1, 
/c = 0, 1, . . . , n2: — 1, 

(Al) 

where 

^ii,m.ini yii,inini ^iZ,min are the minimum values 
of the coordinates in each direction on level il and 
Axii, Ayu^ Azii are the mesh sizes in each direction on 
level il. 

The mesh size between two consecutive levels differs 
by a factor of two so that 



Ax, 



Axi, 



.1+1 



(A2) 



and similarly in the y and z directions. 

In order for the grids to be properly nested we demand 
that there exists an z G [0, — 1] such that 

3 

Xii^i = Xi^+i,min - -Axii^i, il = 0,...nl-2 (A3) 

and similarly in the y and z directions. This condition 
ensures that two consecutive levels share a common in- 
terface. 

We now borrow FMR terminology to name two types 
of cells that exist on our grid. These are the split cells 
and the leaf cells or leaves. A split cell is one within 
which there exist higher level cells and a leaf cell is one 
within which there are no higher level cells. The total 
number of cells A/tot on our grid is 
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A^tot = nx • ny • nz • nl, 



(A4) 



and a straightforward calculation shows that the number 
of leaves is 



Meaf = nx - ny ■ nz 



7n/ + l 



(A5) 



When nl 
expected. 



1, A/'ieaf = A/tot, i-e., all cells are leaves, as 
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We distinguish between these two types of cehs be- 
cause our solutions are defined only on leaves. This may 
be more cumbersome to implement, but has two main 
advantages. 

First, there is no ambiguity as to how one should inter- 
polate values of matter sources from fine levels on coarse 
levels in order to correctly calculate the gravitational 
fields. To be more specific, let us assume that we have 
one coarse cell which is split into 8 cells and that we know 
the values of the density on the fine cells. In Newtonian 
physics, to ensure that the gravitational fields are com- 
puted correctly (at least far away), all we have to do is 
set the cell averaged density on the coarse level such that 
the total mass in the coarse cell is the same as the total 
mass in the enclosed fine cells. In general relativity the 
definition of gravitational mass depends not only on the 
density, but also on the gravitational fields. Hence, there 
is no straightforward way to set the density on coarse 
cells in GR. The ambiguity is immediately lifted, if one 
defines all fields only on the finest cells. 

The second advantage of using only leaves is that the 
memory requirements are minimized and the calculations 
are carried out faster because 



where 



Meaf 7nl + 1 



Snl 



< 1, 



(A6) 



where the last inequality holds because nl > 1. 

For a general second-order nonlinear elliptic equation 
of the form 



V^u = f{u)x, 



(A7) 



where f{u) is a nonlinear function of the variable u and 
X a known scalar independent of li, our code employs a 
standard second-order finite difference scheme 







^'^il,ij,k 










K Uil,i,j-l,k 


^"^11,1 J, k 












'^'^il,i,j,k 





f{Uii^ij^k)Xil,iJ,k' 

(A8) 

The finite difference stencil changes only across grid- 
level boundaries where we perform first order interpola- 
tion. 

To address the nonlinearity of Eq. (|A7p we perform 
Newton- Raphson iterations as follows. Let us assume 
that Un is a guess at step n. We first calculate the residual 
Rn from Eq. (jATl) 



Rn = V'^Un - f(Un)X^ 



(A9) 



and then solve the linearized equation for the correction 
5un on Un 



dm 

du 



(All) 



Once a solution to ([AlOp is found, we correct Un as 

Un+l = SUn, (^12) 

and iterate until this procedure converges and a solution 
to Eq. ()A7p is obtained. 

We solve the equations using the PETSc linear solver 
Krylov Space (KSP) methods. KSP methods are matrix 
methods and hence we have to set up the matrix of the 
linear system. 

To do this we define a global index that counts all cells 
(both leaves and split cells) on our grid as 



il -\- nl{i -\- nx ■ j -\- nx ■ ny ■ k) ^ 



(A13) 



In this way, every leaf cell corresponds to a unique index 
/. However, / takes values 0, 1, ... , A/'tot — 1, but there 
are A/'ieaf leaves on the grid with A/'ieaf < A/'tot- Hence / 
cannot be used to count leaves, unless nl = 1. For this 
reason, we define another index ic, which counts only 
the leaves on our grid, as well as two mappings; from I 
to ic, ic{I) and from ic to I, I{ic). Since for every cell on 
our grid we can find / from Eq. ()A13p we set up these 
mappings by defining two arrays ic_of_I, I_of_ic, of length 
A^tot and Meaf, respectively. Looping over /c, we 

store ic in the array ic_of_I assigning a value of —1 for 
split cells, whereas we store / in the array I_of_ic. The 
index / is used when we need the index ic of a neighbor 
leaf cell in order to calculate derivatives or enter matrix 
elements. 

For example, let us assume that we are at a leaf cell of 
index ic which is not near a grid-level boundary, and we 
want to enter the element of matrix A that corresponds 
to the right-x neighbor of this cell (where A represents 
the Laplacian) . From Eq. ()A8[) this matrix element must 
be 1/Ax|. If ic represents the ic-th row of A we must 
find which column of A the neighbor corresponds to. We 
find this as follows. 

First, using the mapping from ic to I, we find the index 
/ = I{ic) of the leaf. Next by use of Eq. ()A13p we de- 
termine il^i^j^k that correspond to / using the following 
sequence of operations 

k = int{I/nl • nx ■ ny) 

11 = I — nl ■ nx ■ ny ■ k 
j = int{Ii/nl ■ nx) 

12 = h — nl ■ nx ■ j 
i = int{l2/nl • nx) 

il = I2 — nl ■ z. 



(A14) 



V^<5w„ = f'{un)x5u„ - Rn, 



(AlO) 



where int means the integer part of the division. 

In the next step the global index (/pi) of the right x- 
neighbor is found, as /pi = il-\-nl{i-\-l-\-nX'j-\-nX'ny-k). 
Finally, using the mapping from I to ic we find the leaf (or 
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desired column) number ic^i = ic(/pi). Knowing the col- 
umn number of the neighbor, it is straightforward to as- 
sign Aic^icpi = We use the same approach to set 
up all the matrix elements of the linear systems and cal- 
culate derivatives. The algorithm becomes slightly more 
complicated when the neighbor cell is a fictitious cell that 
resides on a different level. In such cases we perform first 
order interpolation and use the same method as outlined 
above to find which matrix elements must be filled with 
non-zero values. 



Appendix B: Validation of HRSC method for a 
piecewise polytropic EOS 

In this appendix we analyze the effect our numerical 
schemes have on solutions obtained with our adopted 
non-smooth EOS (fT5]) and a smooth counterpart of this 
EOS. We show that there is no essential difference. This 
result is expected because an algorithm with finite resolu- 
tion cannot distinguish a smooth EOS from a non-smooth 
EOS, if the smoothing operation is performed below the 
resolution level of the computations. 

For simplicity we consider a non-smooth EOS with two 
branches as follows 



l+l/m ^ 



l + l/ns 



(Bl) 



and perform a smoothing operation over a density inter- 
val — e), + e)] as follows 



l+l/ni ^ /i \ 



/(Po), 



pi(l-e)<po^Pi(l + e) , (B2) 



l+l/n2 . /i I \ 

I i^2Po , Po > Pi(l + e). 



where /(po) is a smooth spline fit such that the EOS is 
continuous and has continuous first or second derivative, 
depending on the order of the spline. Our particular 
choice for the smoothing function is either a cubic spline 
or a quintic spline. In the former case the EOS becomes 
C^, while in the latter case the EOS becomes C^. 

We chose e to be sufficiently small so that the smoothed 
EOS mimics as closely as possible EOS (|Bip . but large 
enough to avoid round off errors due to very large gra- 
dients. For the cubic spline we set e = 10~^, while 
for the quintic spline we set e = 10~^. In all our nu- 
merical tests we choose /ci, A:2, ni, n2, pi to correspond to 
/c2, /c3, 712, na, p2 of the 10:1 EOS (see Table |I|), respec- 
tively. In Fig. [13] we show a plot of EOSs (|Bl]) and (IB2]) . 
where /(po) is a cubic spline. 

With the smooth EOS ()B2p at our disposal we set up 
several ID Riemman (or shock tube) problems in a spa- 
tial domain of length L = 4.2km and resolutions 210, 
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FIG. 13: Pressure vs rest-mass density for EOSs ()B1|) (black) 
and ^^^-^^ "^^^^ - 

IQ-^km" 

non-differentiable, and shows that the cubic spline fit smooths 
out the discontinuity. 



(red). Here Po = lO^'^km"^ and pnuc = 1.48494 x 
^. The inset zooms in the region, where EOS ()B1|) is 



420, 840 and 1680 grid points. We set Fth = 1.66 and 
use Eq. (j40]) with Pcoid given by either EOS (jBl| or EOS 
(|B2l) . 

We have explored the parameter space of initial data 
{Pl^ Pl^Ul)^ {Pr^ Pri^r) the left and right states, 
and in all cases we found that the solutions obtained 
with EOS (|Bip almost overlap those obtained with the 
smooth EOS (|B2p . These results hold for both piece- 
wise parabolic (PPM) and monotonized central (MC) re- 
construction, regardless of resolution and the spline fit 
choice. Furtermore, we verified that all our simulations 
with the smooth EOS ()B2p had high enough resolution 
so that data points would sample the smoothing interval 
[pi(l — e), pi(l + e)] every few timesteps. 

In Fig. [m we plot a snapshot of the pressure profile 
and do a convergence test for one of the cases we sim- 



ulated with {pR = 10 Pr = Pcoid(pi?), 



= 0), 



(pL = 5 X 10-"^, Pl = IOPr, ul = 0). The figure shows 
that all solutions overlap (left panel) and converge at 
the expected order (right panel) to the very high reso- 
lution solution obtained with PPM in conjunction with 
the smooth EOS. 

The solutions obtained with the smooth EOS ()B2p with 
quintic spline smoothing, which results in a and a 
convex EOS, also overlap with those of the non-smooth 
EOS solutions even though the smoothing interval we 
chose was much larger than in the cubic spline case, and 
hence the data points would sample it more frequently. 
Note also that the coarsest resolution used in the shock 
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FIG. 14: Left: Snapshot of pressure profile at t = 2.56 km, corresponding to a tiiird of a sound speed crossing timescale 
across the domain. Right: Convergence plot (at t — 2.56 km) using as reference solution the very high resolution solution 
obtained in conjunction with the smooth EOS ()B2|) with a cubic spline smoothing function. The labels in the plots denote 
the resolution (LR, MR, HR, or VHR) and the reconstruction method (PPM or MC as subscripts). The resolutions used are: 
LR = 210 , MR = 420, HR = 840, VHR = 1680 grid points. PPM stands for the piecewise parabohc reconstruction, and MC 
stands for the monotonized central reconstruction. The plots demonstrate that all solutions overlap (left panel), regardless of 
the reconstruction method and the EOS used (smooth or non-smooth), and first-order convergence to the VHR run with the 
smooth EOS, as expected. Here Pq — 10~^km~^. 



tube tests is at least 20 times higher than the resolu- the correct solution and indicate that no unphysical solu- 
tion used in our WDNS simulations. Therefore, all these tions are present in our simulations of the WDNS head-on 
results demonstrate that our numerical methods capture collisions. 
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